
#ifndef HJM_HPP
#define HJM_HPP

#include "Matrix.hpp"
#include "Normal.hpp"
#include "ThreadedObserver.hpp"

#include <boost/thread.hpp>
#include <boost/thread/mutex.hpp>
#include <boost/random.hpp>

using namespace std;
using namespace boost;

// Mersenne twister
typedef boost::mt19937 base_generator_type;

class HJM : public ThreadedObserverSubject {
  vector <double> m_X; // maturities
  vector <double> m_M; // drift
  Vector <double> m_F0; // initial rates
  Matrix <double> m_EigenVec; // volatility component
  vector <double> m_dX;
  //
  // simulation results
  //
  vector < Vector <double> > m_TS; // time series
  vector < Matrix <double> > m_F; // simulated rates
  vector < Vector <double> > m_SRate; // short rate
  vector < Vector <double> > m_Disc; // discount
  //
  // Mutex for threads reporting results.
  //
  mutex m_mutex;
  //
  // Simulation params
  //
  int m_nPaths; // # paths to simulate
  int m_nSteps; // # steps in simulation
  int m_nMat; // # maturities
  int m_nThreads; // # threads
  double m_Tf; // final simulation time
  double m_dt; // time step
public:
  HJM ( const vector<double>& X, const vector<double>& M, 
	const Matrix<double>& p_EigenVec, const Vector<double>& F0, 
	double p_Tf, int p_nSteps, int p_nPaths, int p_nThreads );

  virtual ~HJM();
  void simulate (int ith, ulong p_seed);
  //const Matrix<double>& Rates() const { return m_F; }
  //const vector<double>& TimeSeries() const { return m_TS; }

  double Caplet ( int ith, double strike, double tenor, double Tf );
  double Floorlet ( int ith, double strike, double tenor, double Tf );
  // discounted sum of series of caplets
  //double cap ( double strike, double tenor, double Tf );

  inline double ZCB ( int ith, double Tf )
  { 
    return Discount(ith,Tf);
  }
  // calculate the discount rate
  inline double Discount ( int ith, double Tf )
  {
    vector <double>& t_TS = m_TS[ith];

    int nT = t_TS.size();
    if ( t_TS[nT-1] < Tf )
      throw (range_error("HJM::Discount: Tf out of range."));

    // Tf - find the closest
    double eps = 1.0e-4;
    int kf = 0;
    while ( t_TS[kf] < Tf && abs(Tf - t_TS[kf]) > eps )
      ++kf;

    return exp ( -m_dt * m_Disc[ith][kf] );
  }
  // Libor - mean rate over an interval
  double Libor ( int ith, double Tmat, double Tstl );
  double Call ( int ith, double strike, double Tf );
  double Swap ( int ith, double period, double Tf );
};

#endif
